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Abstract 

We derive a new completely integrable dispersive shallow water equation that is 
biHamiltonian and thus possesses an infinite number of conservation laws in invo- 
lution. The equation is obtained by using an asymptotic expansion directly in the 
Hamiltonian for Euler's equations in the shallow water regime. The soliton solution 
for this equation has a limiting form that has a discontinuity in the first derivative 
at its peak. 
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Completely integrable nonlinear partial differential equations arise at various levels 
of approximation in shallow water theory. Such equations possess soliton solutions - 
coherent (spatially localized) structures that interact nonhnearly among themselves then 
re-emerge, retaining their identity and showing particle-like scattering behavior. In this 
paper, we use Hamiltonian methods to derive a new completely integrable dispersive 
shallow water equation. 



where u is the fluid velocity in the x direction (or equivalently the height of the water's 
free surface above a flat bottom), «; is a constant related to the critical shallow water 
wave speed, and subscripts denote partial derivatives. This equation retains higher order 
terms (the right-hand side) in a small amplitude expansion of incompressible Euler's 
equations for unidirectional motion of waves at the free surface under the influence of 
gravity. Dropping these terms leads to the Benjamin- Bona- Mahoney (BBM) equation, 
or at the same order, the Korteweg-de Vries (KdV) equation. Our extension of the 
BBM equation possesses soliton solutions whose limiting form as k ^ have peaks 
where flrst derivatives are discontinuous. These "peakons" dominate the solution of the 
initial value problem for this equation with k = 0. The way a smooth initial condition 
breaks up into a train of peakons is by developing a verticality at each inflection point 
with negative slope, from which a derivative discontinuity emerges. Remarkably, the 
multisoliton solution is obtained by simply superimposing the single peakon solutions 
and solving for the evolution of their amplitudes and the positions of their peaks as a 
completely integrable finite dimensional Hamiltonian system. 

Our equation is biHamiltonian, i.e., it can be expressed in Hamiltonian form in two 
different ways. The ratio of its two (compatible) Hamiltonian operators is a recursion 



Ut + KUx — Uxxt + 3mM; 



'a; 



'XXX ) 
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operator that produces an infinite sequence of conservation laws. Tliis biHamiltonian 
property is used to recast our equation as a compatibility condition for a linear isospec- 
tral problem, so that the initial value problem may be solved by the inverse scattering 
transform (1ST) method. 

The unidirectional model. Consider Euler's equations for an inviscid incompressible 
fluid of uniform density with one horizontal velocity component u in the x direction, and 
w in the vertical {z) direction. The fluid is acted on by the acceleration of gravity, g, 
and is moving in a horizontally infinite domain with an upper free surface at z = C{x,t) 
and flat bottom at z = —ho. Substituting the solution form motivated by shallow water 
asymptoticsQ , u = u{x,t), w = —{z + ho)ux, into the conserved energy (kinetic + 
potential) for Euler's equations, and explicitly performing the 2;-integration leads to the 
energy Hgn = \ J-^ dx rju^ + \rful. + gij] — , where ?7 = C + is the height of 
the water above the bottom. Substituting the same solution form above into Euler's 
equations and integrating over the vertical coordinate leads to the Green-Naghdi (GN) 
equation J. The GN equation, conserve the energy In tact, they ate expressible 

in Hamiltonian form© as 



rrtt\ _ ( dm + md rjd \ f 5Hgn/^^ 
Vt dr] oil 6Hgn/Sv 



(2) 



where the momentum density m is defined by m = 6Hgn/Su. The GN equations do not 
necessarily refer to a thin-domain expansion in a small parameter e that measures the 
ratio of depth to wavelength. In such an expansion the kinetic energy of vertical motion 
(~ Tj^ul) in Hgn would be O(e^). Shallow water theory makes a further small- amplitude 
assumption, in the form r] = ho + 0{a), a « 1, and balances a = O(e^). In contrast, 
the Hamiltonian Hgn retains nondominant terms (e.g., that would be higher order in 
such an expansion. Starting from the GN equations, further small-amplitude asymptotics 
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and restriction to unidirectional propagation in a frame moving near the critical wave 
speed Co = \/gfk>, leads to the KdV equation 0, + cou^ + 3/2 uu^ + 1/6 Cq/iqW^xx = 0, 
or, with the same order of accuracy in the thin-domain expansion, the BBM equationi, 
Ut + CqUx + 3/2 uUx — 1/6 h^Uxxt = 0. Instead of making asymptotic expansions in 
the equations of motion, as in the derivations of the KdV and BBM equations, our 
approach in deriving (|1]) is to make a unidirectional approximation by relating m to rj 
in the GN system and preserving the momentum part of its Hamiltonian structure (0). 



For this purpose, we will set r] = hQym/{hoCo), and since rj ^ ho as \x\ oo the 
boundary conditions on m will be assumed to be m ^ /iqCq as |a;| —>■ oo. The functional 
= V^n dx is the Casimir for the Hamiltonian operator {md + dm) and so we will 
refer to this invariant manifold as the Casimir manifold for (^. Next, we scale m — > aw in 
the Hamiltonian Hon, look for m in the form m = /iqCq + ami + (y'^m2 + a'^m^ + . . . and 
expand Hgn accordingly. With this scaling and expansion, defining m as the variational 
derivative of tire Hamiltonian with respect to and balancing at order O(c^) givesi 
mi = 2{hQU — h^Uxx/'i)- The Hamiltonian may then be rewritten as i^cAf = Hi£, + 
0{a^), where Hid = «^/4 miu dx + a/2 miCo dx, and the factor 1/2 arises 
from restricting to a submanifold.01 

The 0{a) equation of motion for m on the Casimir manifold is therefore mt = —{md+ 
dm)5Hi£)/ 5m = —a/2{md + dm)u — Cq/2 m^, or, in terms of u, 

1 3 111 

Ut - -hlu^xt + CoMa; + -auu^ - -hlcQU^xx = -ahlu^u^^ + -ahluu^^^. (3) 
3 2 6 3 6 

Dropping the right-hand side of this equation gives BBM or KdV, modulo replacing 
Uxxt by — Co'Uxxx-'S Thus (H) can be seen as a BBM equation extended by retaining 
higher order terms (selected by the Hamiltonian approach) in an asymptotic expansion 
in terms of the small-amplitude parameter a. The restriction to the Casimir manifold 



Peakons, Camassa & Holm — (February 7, 2008) 



5 



is equivalent at order 0{a) to the unidirectionality assumption ( = ^Jho/g u + 0{a) in 
the usual derivations of the KdV and BBM models from the Boussinesq system'S'S. In 
fact, C = \JhQ/ g [u — hyQux^^ + 0{a), and in a thin-domain approximation the double 
derivative term in this expression would acquire a factor e^. 

Rescaling (^), dropping a, and going to a frame of reference moving with speed 
n = Co/2 reduces the equation to the standard form (|I|). Notice that (|T]), like BBM, is 
not Galilean invariant, i.e., not invariant under u ^ u + k, t ^ t, x ^ x + nt. Thus, 
equation (|^) is best seen as a member of a family of equations parameterized by the speed 
K of the Galilean frame. 

Using the identity (1 — 9^)e~l^l = 25{x) and setting }C[v] = dyex-p{—\x — y\)v{y), 
expresses equation (|T]) in nonlocal form as Ut + uu^ + i^lC[uy\ = —IC[uUy + 1/2 UyUyy]. 
Dropping the quadratic terms on the right-hand side of this equation gives the one studied 
by Fornberg and WhithamS . Fornberg and Whitham show that traveling wave solutions 
of this truncated equation have a peaked limiting form. Moreover, nonsymmetric initial 
data with two inflection points in their case can develop a vertical slope in finite time. 

In a later paper we will discuss the parameterized family (|T]). The present paper 
focuses on the limiting case k = 0, 

Ut - U^xt = -SuUx + 2UxUxx + UUxxx , (4) 

where u is defined on the real line with vanishing boundary conditions at infinity and such 
that the Hamiltonian Hi = ^ + "^D bounded. As with (|^), Hi generates 

the flow (^) through m = u — u^xi TTtt = —{md + dm)5Hi/ 5m. 

Steepening at inflection points. Consider an initial condition that has an inflection 
point at X = X, to the right of its maximum, and decays to zero in each direction suffi- 
ciently rapidly for Hi to be finite. Define the the time dependent slope at the infiection 
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point as s{t) = Ux{x{t),t). Then the nonlocal form of (^) (with k = 0) and standard 
Sobolev estimates yield a differential inequality for s, ds/dt < — + Hi. Hence, 
the slope becomes vertical in finite time, provided it is initially sufficiently negative. If 
the initial condition is antisymmetric, then the inflection point at m = is fixed and 
dx/dt = 0, due to the symmetry {u,x) {—u,—x) enjoyed by (H). In this case, no 
matter how small |s(0)|, verticality develops in finite time. This steepening property im- 
plies that traveling wave solutions of cannot have the usual bell shape since inflection 
points may not be stationary in time. In fact the traveling wave solution is given by 
u{x,t) = cexp(— |x — ct\). This solution travels with speed c and has a corner (that is, a 
finite jump in its derivative) at its peak of height c.0 

A^-soIiton solution. Motivated by the form of the traveling wave solution, we make the 
following solution ansatz for interacting peaked solutions, u{x, t) = Y^iLi Pi{^) exp(— |x — qi{t)\). 
Substituting this into equation (^) yields evolution equations for and pj, that are 
Hamilton's canonical equations, with Hamiltonian Ha given by substituting the solution 
Ansatz above into the integral of motion Hi , yielding Ha = 1/2 Yl,fj=i PiPj exp (— \qi — qj\)- 
Hamiltonians of this form describe geodesic motion. The peak position qi{t) is governed 
by geodesic motion of a particle on an A^- dimensional surface with inverse metric tensor 
9^^{q) = 6xp(— |gj — 9 G R^. The metric tensor is singular whenever = qj. 
Two-soliton Dynamics. Consider the scattering of two solitons that are initially well 
separated, and have speeds Ci and C2, with ci > C2 and ci > 0, so that they colhde. The 
Hamiltonian system governing this collision possesses two constants of motion, Hq = 
Pi + P2 = Ci + C2 and H^ = {cf + cl)/2. Notice that if the peaks were to overlap, 
thereby producing qi — q2 = during a collision, there would be a contradiction 2Ha = 
(ci + €2)^ = Ci + C2, unless p were to diverge when the overlap occurred. 

The solution of Hamilton's canonical equations for Hamiltonian H^ when A^ = 2 is 
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given by 

4(ci - C2)27e(^i-"2)t 



gi - 92 = - log 



Pl-P2 = ±(Cl-C2) 7 7-. (5) 

and the conservation law for pi+p2- Here 7 is a constant specifying the initial separation 
of the peaks, and ci and C2 are the asymptotic t —>■ ±00 values of their speeds, or 
amplitudes. The divergence of pi and p2 in equation (|^) associated with soliton overlap 
can only occur when ci and C2 have opposite signs. That is, only "head-on" collisions can 
lead to overlapping peaks (see Fig. 1, available from authors, for the "soliton- antisoliton" 
case ci = — C2 = c). 

The two soliton solution determines the "phase shifts," i.e., the shifts in the 
asymptotic position for t —>■ 00, that the solitons experience after interaction. Ast +00 
the solitons re-emerge unscathed, the faster (and larger) soliton ahead of the slower (and 
smaller) one. Defining the phase shift for the faster soliton to be Aqf = g2(+oo)— gi(— 00), 
and for the slower soliton, Aqg = gi(+oo) — q2{—oo), leads to Aqj = log(cf/(ci — C2Y), 
and Aqs = log((ci — C2Y/C2). These formulae show that when Ci/c2 > 2 both solitons 
experience a forward shift. For 1 < Ci/c2 < 2 the faster soliton is shifted forward while 
the slower soliton is shifted backward. When Ci/c2 = 2 no shift occurs for the slower 
soliton. 

BiHamiltonian structure. Equation (^) follows, as well, from an action principle 
expressed in terms of a velocity potential. This action principle leads to an additional 
conserved quantity, H2 = \ I-^{u^ + uu^dx, and another Hamiltonian operator, d — 
d^. Our equation (^ then can be written in Hamiltonian form in two different ways, 
rrit = —{d — d^)5H2/ 5m = —{md + dm)5Hi/ 5m. The two Hamiltonian operators Bi = 
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d — d^, and B2 = dm + md form a Hamiltonian pair. That is, their sum is still a 



Hamiltonian operatorl^JJ. Equation (^) is thus biHamiltonian and has an infinite number 
of conservation laws recursively related to each other by BiSHn/Sm = B26Hn-i/ 5m = 
—m["'^^\ n = 0, ±1, ±2, . . .. Starting from Hi or H2 this relation generates an infinite 
sequence of conservation laws including, e.g., Hq = m dx, H^i = j^!^ \/m dx = 
C, H_2 = 1 1-^ m^/4m^/^ — 2/y/m dx, etc. Correspondingly, the recursion operator 
TZ = i?2-Br^ generates a hierarchy of commuting flows, defined by ml"^^^ = Kn+i[m] = 
TZKn[m\, n = 0, ±1,±2,.... The first few flows in the hierarchy are mf^^ = —{d — 
^^){2^/m)^^ , m[^^ = 0, m[^^ = —m^, and mp"* = —{md + dm)u. The last of these is our 



equation and the first is an extension of the integrable Dym equationOl . It turns 
out that all the flows in this hierarchy are isospectral and thus completely integrable. 
The isospectral problem. In order to find the isospectral problem for our equation, we 



follow Gel'fand and Dorfmannli3 in considering the skew symmetric spectral problem, 
{\Bi — -82)0 = 0. A class of solutions of this problem are related hy (f) = ip"^ to the 
solutions ip oi a, second order symmetric spectral problem. By imposing isospectrality. 
At = 0, our equation (|l]) follows from the compatibility condition ipxxt = '4'txx of the 
system for ^/^(x, t). 



'^x 



1 m{x, t) -\- n 



^, ipt = -{\ + u)iJx + ^Uxij. (6) 



4 2A 

This is the isospectral problem we seek. The system (|^) provides a means of solving the 
initial value problem for (|1|) by the purely linear 1ST technique&3. For instance, if the 
boundary conditions on m are taken to be zero at x = ±00 (sufficiently fast),0 then the 
spectral problem when k = has a purely discrete spectrum since ip{x) exp(±x/2) 
as |x| 00, i.e., eigenfunctions always decay exponentially at infinity. If, e.g., the initial 
condition u{x,0) is chosen such that u{x,0) = A{7t/2 — 2 sinhx arctan (e^) — 1), so 
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that m{x,0) = A sech^(x), for an arbitrary constant A, then it is easy to showC-3 that 
the eigenvalues A for (||) are given by A„ = 2A/[{2n + l){2n + 3)], = 0, 1, 2, . . .. This 
formula shows explicitly that A = is an accumulation point for the discrete spectrum 
and the eigenvalues converge to it as n oo, a fact that holds in general for 

any initial condition decaying exponentially fast at infinity. Equations also imply 
that the A^-soliton mechanical system with Hamiltonian Ha is completely integ rableEl. 
When K 0, i.e., for an equation in the family ([l|), the limiting behavior of ip becomes 



ipi^x) — > exp (^±x^l/4: — k/2X^ as |a;| — > oo, and so continuous spectrum develops out of 
the origin in the interval < A < 2k. Also, for k 7^ the soliton solution of (|l|) becomes 
C°° and there is no derivative discontinuity at its peak. The peculiar feature of the 
disappearance of continuous spectrum in the limit k can be traced to the constant 
1/4 in the spectral problem (P), which in turn is generated by the first derivative operator 
in Bi. 

Numerical .imulation^H confirm the a„aly=i= discussed here and demonstrate the 

robustness of the peaked soliton solutions. These simulations clearly illustrate the inflec- 
tion point mechanism by which a localized (positive) initial condition breaks up into a 
height-ordered train of peaked solitons moving to the right, with the tallest ones ahead. 

For their helpful and sometimes challenging remarks during the course of this work, 
we thank M. Ablowitz, I. Gabitov, I. M. Gel'fand, J. M. Hyman, B. Kupershmidt, P. 
Lax, C. D. Levermore, S. V. Manakov, L. Margolin, P. Olver, T. Ratiu, H. Segur and T. 
Y. Wu. This work is partially supported by the U.S. Department of Energy CHAMMP 
program. 



Figure caption. 
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Fig. 1. The soliton-antisoliton solution u reconstructed from equation (^) is u{x, t) = c[exp(— |x— 
1/2 q(t)\) — exp(— |rE + 1/2 g(t)|)]/tanh(ct). This solution displays the steepening behavior dis- 
cussed in the text. The slope becomes vertical and the amplitude of the solution becomes 
(everywhere) zero right at the moment of overlap. At later times the peaks redevelop and 
depart again according to the symmetry {u,t) — > {—u, —t). 
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